%varying lamdac
clear all
clc
muc=0.033;
C=10:1:40;
Lambdac= 5 ;
%Lambdac=3;

Lambdav=8;
Waitings=zeros(length(C),1);
R=zeros(length(C),1);
exitflags=7*ones(length(C),1);
options = optimset('MaxFunEvals',40000);
for j=1:length(C)
n=ceil(((Lambdav/muc) -1)/C(j));
% Pv=Gaussiandist(n);
Pv=decreasingdist(n);
Pc=Gaussiandist(n);
var_num = ((n+1)*(n+2))/2; 

lb=zeros(var_num,1);
ub=ones(var_num,1);
ub(var_num)=Inf;
A1=zeros(1,var_num);
for i=1:n
A1(1,i)=-Lambdav*Pv(1,i);
end
b1=ones(1,1)*muc*n*C(j) - Lambdav;
A2=zeros(1,var_num);
A2(1,1)=Lambdav*Pv(1,1);
b2=ones(1,1)*muc;
Alin=[A1;A2];
blin=[b1;b2];

Aeq=zeros(n,var_num);
beq=ones(n,1);
begins=n;
ends=n;
for i=1:n
    begins=ends+1;
   ends=ends+i;
    for k=begins:ends
        Aeq(i,k)=1;
    end
   
end

b=zeros(n,1);

for i=1:n
    b(i)= Lambdac*Pc(i) ;
end
ub(var_num)=(Lambdav-Lambdac)/n;
cond = @(x)parameterfun(x,n,var_num,Pv,Lambdav,b);
obj=@(x)x(var_num)*(-1);
x0=zeros(var_num,1);
%x0 = ones(var_num-1,1);
%x0(var_num)= (Lambdav-Lambdac(j))/n;
x_best=zeros(var_num,1);
    for q=1:20
        [x,fval,exitflag,output,lambda] = fmincon(obj,x0,Alin,blin,Aeq,beq,lb,ub,cond,options);
        %x(var_num)
        if x(var_num)> x_best(var_num)
            x_best=x;
            exitflag_best=exitflag;
            lambda_best=lambda;
            %best_result=[x,fval,exitflag,output,lambda];
            x0=x;
        end
    end
exitflags(j)=exitflag_best;
R(j)= (1/x_best(var_num));
%[cout,ceqout] = parameterfun(x,n,Pv,Lambdav,b);
%Waitings(j)=max(1/(-cout(1:n)+0.1));
end


plot(C,R)
hold on